function [] = trackingSystem()

    theta = 0; % actual angle of the radar antenna ...
    theta_R = []; % angle between target and the radar antenna ...
    
    %K = [10, 5, 1, 0.5, 0.3, 0.1, 0.01, 0.001, 0.0001]; % constant values of the torque controller for the antenna ...
    %K = [0.1, 0.01, 0.005]; % constant values for the integrative controller ...
    K = [-0.1, -0.03, -0.01, -0.001]; % negative constant values to demonstrate instable systems (derivative contr.) ...
    ki = 0; % temporary variable for the constant values of the torque controller ...

    %% constant parameters for the antenna:
    IM = 0.004; % inertia moment of the antenna in [kg m^2] ...
    b = 0.02; % viscosity coefficient [kg m^2 s^-1] ...
    
    %% nested functions:
    %% antenna models:
    
    % proportional model of the antenna:
    function x_prim = propAntennaModel(ta, ya)%, ki)
        x1 = ya(1);
        x2 = ya(2);
        theta_r = bearTarget(ta);
        
        % coefficient matrix A ...
        coeffMat_A = [0, 1, 0; -ki/IM, -b/IM, (ki*theta_r)/IM];
        % vector with unknown variables ...
        vec_x = [x1; x2; 1];
        % create the model of the system by applying the formula
        % dx/dt = A*x + b:
        %(note: here is b embedded in matrix A and vec_x.)
        x_prim = coeffMat_A * vec_x;
    end

    % integrative model of the antenna:
    function x_prim = intAntennaModel(ta, ya)
        x1 = ya(1);
        x2 = ya(2);
        x3 = ya(3);
            
        % coefficient matrix A ...
        coeffMat_A = [0, 1, 0; 0, 0, 1; -ki/IM, 0, -b/IM];
        % vector with unknown variables ...
        vec_x = [x1; x2; x3];
        % vector with constant values ...
        vec_b = [0; 0; (ki*0.005*ta^2)/IM]; % 0.01/2 = 0.005
        % calculate the equation system (dx/dt = A*x + b):
        x_prim = coeffMat_A * vec_x + vec_b;
        
        % correction: swap the first two rows of the vector (why?) ...
%         x1_prim = x_prim(1,:);
%         x2_prim = x_prim(2,:);
%        
%         x_prim(1,:) = x2_prim;
%         x_prim(2,:) = x1_prim;
    end
    
    % derivative model of the antenna:
    function x_prim = derivAntennaModel(ta, ya)
        % note: we don't need here the time values "ta" because of dx/dt ...
        x1 = ya(1);
        x2 = ya(2);
        
        % coefficient matrix A ...
        coeffMat_A = [0, 1 ; 0, -(ki+b)/IM];
        % vector with unknown variables ...
        vec_x = [x1; x2];
        % vector with constant values ...
        vec_b = [0; 0.01*ki/IM]; 
        % calculate the equation system (dx/dt = A*x + b):
        x_prim = coeffMat_A * vec_x + vec_b;        
    end

    %% bearing:

    % method to bear the target (calculate the position = angle):
    function theta_r = bearTarget(t)
        theta_r = 0.01*t;
    end

    %% error methods:
    
    % calculates the signal error e(t):
    function e = signalError(theta, theta_R)
        %e = abs(theta_R - theta);
        e = theta_R - theta;
    end

    % relative error of the signals:
    function re = relError(theta, theta_R)
        re = signalError(theta, theta_R)./theta_R;
    end


    %% equation for tracking the target:
    %    IM*Theta'' = -b*Theta' + u(t)
    %    u(t) = K*e(t), e(t) = difference between "theta_R" and "theta"

    %% simulation:
    
    % initial condition:
    y0 = [0; 0]; % proportional and derivative model ...
    %y0 = [0; 0; 0]; % integrative model ...
    
    t_i = 0;        % initial value
    t_f = 300;      % final value
    %iNSteps = 500; % > 450 steps are need to get the same result as with "ode45" ...
    
    % file handler to the model function ...
    %fhandle = @propAntennaModel;
    %fhandle = @intAntennaModel;
    fhandle = @derivAntennaModel;

    clf; % clear the frame of the plot ...
    % create a color list (grey colors for each line) ...
    cGray = [0.4, 0.4, 0.4];
    grayColor = [ [0.1, 0.1, 0.1]; [0.2, 0.2, 0.2]; cGray; ...
                  [0.5, 0.5, 0.5]; [0.5, 0.5, 0.5]; ...
                                
                  % colors for the grafic of question a.) and e.) - inst.
                  % systems ...
                  [0.5, 0.5, 0.5]; [0.5, 0.5, 0.5]; ...
                  [0.5, 0.5, 0.5]; cGray ];
                  
                  % colors for the grafic of questions b.), d.) and e.) -
                  % stable systems ...
                  %[0.6, 0.6, 0.6]; [0.6, 0.6, 0.6]; ...
                  %[0.7, 0.7, 0.7]; [0.7, 0.7, 0.7] ];
    lineType = {'--', '-.', ':', '--', '-.', ':', '--', '-.', ':'};
    
    %bearVals = zeros(length(K), (iNSteps + 1));
    
    %% solving the ODE and plot the results:
    
    for i = 1:length(K) % loop for the other questions ...
    %for i = 7:length(K) % loop for question a.)
    %for i = 2:length(K) % loop for question e.) - stable systems 
        ki = K(i);
        
        %[t, y] = solveODE_RK4(fhandle, t_i, t_f, y0, ki, iNSteps); % buggy ...
        %[t, y] = solveODE_RK4(fhandle, t_i, t_f, y0, iNSteps);
        %bearVals(i,:) = y(1,:);
    
        [t, y] = ode45(fhandle, [t_i, t_f], y0); % RK-method (supported by matlab) ...

        % angle curves for question a.)
        % plot of different curves of the angle values approximated with
        % different constant values ...
        plot(t, y(:, 1), lineType{i}, 'Color', grayColor(i,:), 'LineWidth', 1);
        hold on;
        
        % error curves for question c.)
        %relerror = relError(y(:, 1),t.*0.01);
        %plot(t,relerror);
        %hold on;
    end
        
    % calculate the (original) angle values between target and the radar antenna ...
    theta_R = bearTarget(t);
    %disp(theta_R');
    
    % plot the original angel values between the target and the antenna ...
    plot(t, theta_R, '-k', 'LineWidth', 1);
%     axis([0 300 0 3]); % axis config. for question a.) and e.) (only for stable systems)

%     axis([0 1 0 0.01]); % axis config. for question b.)
%     axis([0 0.5 0 0.005]); % axis config. for question b.)

%     axis([0 300 0 100]); % axis config. for quest. d.) - show the oszillations ...
%     axis([0 21 -2 5]); % axis config. for quest. d.) - oszillations more in detail ...
%     axis([0 10 0 0.7]); % axis config. for quest. d.) - oszillations in detail at the begin ...
    
     axis([0 2.5 -0.06 0.06]); % axis config. for quest. e.) - instable system

    axis normal
    xlabel('Tiempo [s]', 'FontName', 'Times', 'FontSize', 11);
    ylabel('Angulo [rad]', 'FontName', 'Times', 'FontSize', 11);
    title('',  'FontName', 'Times', 'FontSize', 12);
        
    % legend for the grafic of question a.)
%     h = legend('k_1 = 0.01', 'k_2 = 0.001', 'k_3 = 0.0001', 'blanco');
    
    % legend for the grafic of question b.)
%     h = legend( 'k_1 = 10', 'k_2 = 5', 'k_3 = 1', 'k_4 = 0.5', ...
%             'k_5 = 0.3', 'k_6 = 0.1', 'k_7 = 0.01', ...
%             'k_8 = 0.001', 'k_9 = 0.0001', 'blanco' );

    % legend for the grafic of question d.) - stable systems ...
%     h = legend('k_1 = 0.1', 'k_2 = 0.01', 'k_3 = 0.005', 'blanco');

    % legend for the grafic of question e.) - stable systems ...
%     h = legend( 'k_1 = 5', 'k_2 = 1', 'k_3 = 0.5', ...
%                 'k_4 = 0.3', 'k_5 = 0.1', 'k_6 = 0.01', ...
%                 'k_7 = 0.001', 'k_8 = 0.0001', 'blanco' );

    % legend for the grafic of question e.) - instable systems ...
   h = legend('k_1 = -0.1', 'k_2 = -0.03', 'k_3 = -0.01', 'k_4 = -0.001', 'blanco');

    set( h, 'Interpreter', 'tex', 'Location', 'NorthWest', ... %'Best', ...
         'FontName', 'Times', 'FontSize', 10 );
    hold off;


%% not necessary ...
%     %% control methods to control the motor torques of the antenna:
%     % simple controller to control the motor torques by using a (deterministic)
%     % constant proportional value "K".
%     function u =propController(theta, theta_R, K)
%         u = K*signalError(theta, theta_R);
%     end
% 
%     % integral controller: u(t) = k_i*int(e(tau), tau, 0, t)
%     %
%     % (--> will getting suck in the simulation ;-) )
%     function u = intController(k_i, t)
%         u = k_i*int(e(tau), 0, t);
%     end
% 
%     % derivative controller: u(t) = k_d*(d e(t)/d t)
%     %
%     % (--> will return definitely the better solutions ...)
%     function u = derivController(k_d, t)
%         u = k_d*diff(e(t), t);
%     end

end

 